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SUMMARY 


Shock-capturing solutions for an axisymmetric supersonic inlet at small 
angles of attack have been obtained. Good overall agreement between the 
shock-capturing solutions and experimental data have been shown except in 
regions of strong viscous effects or boundary-layer removal. Although the . 
results Indicate a strong potential for the use of shock-capturing or finite- 
difference solutions for Internal flows, 'Improvement In the ability to handle 
the reflection of strong shockwaves having downstream Mach numbers near 1 is 
needed. 


INTRODUCTION 


Traditional approaches to the design of axisymmetric , supersonic inlets 
have employed the method of characteristics (refs. 1 to 3). References 2 and 3 
Include engineering modeling of the viscous effects, with reference 3 incor- 
porating the effects of the boundary layer on the Inviscid core flow. For 
two-dimensional flows — Including axisymmetric flows — the method of character- 
istics provides a standard for any other computational technique. However, 
extension of the method of characteristics to three-dimensional flows, particu- 
larly complex Internal flows, presents a very formidable task. 


Recent developments In the use of finite-difference techniques show promise 
for application to three-dimensional Internal flows. References 4 and 5 have 
been concerned with solutions to the flow about shuttle-type vehicles, both 
with diffuse (captured) and discrete bow shockwaves. References 6 and, 7 
describe techniques to fit discrete shockwaves into the computational mesh. 

A comparison of a discrete-shock, finite-difference technique with the method 
of characteristics for planar Internal flows was given in reference 7; agree- 
ment between the two techniques was excellent. An advantage in using. the • 
finite-difference codes referenced above, as opposed to the method of charac- 
teristics, is that embedded shockwaves are captured in the finite-difference 
mesh without requiring any special logic. 


The present paper describes an adaptation of the shock-capturing tech- 
nique of reference 4 to the particular problem of inviscid flow in axisymmetric 
supersonic inlets at small angles of attack. (Many helpful suggestions con- 
cerning the present paper were made by Paul Kutler.) Comparison of the theo- 
retical solutions with experimental data for the ^ = 3.5 inlet described in 
reference 2 is presented. The present shock-capturing solutions represent the 
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first step In a. long range program that will ultimately Incorporate discrete ^ 
shockwaves, boundary- layer effects, mass exchange, and an arbitrary body cross 
section. The computer code for this analysis was developed by Dennis J. Itelne, 
Computer Science Corporation. 

SYMBOLS 

a local speed of sound - ■ . i 



M Mach number 

p pressure normallxed by the f ree-stream stagnation pressure 
q ' total-velocity normalized by the maximum adiabatic velocity 
r radius from axis- of symmetry 

u z-component of velocity 

V r-component of velocity 

w (j>-component of velocity 

Z : axial coordinate from centerbody tip 

a angle of attack 

Y ratio of specific heats 

6 boundary surface angle 
p denslt:y 

^ meridional coordinate 
Subscripts ; 
a annulus 

b centerbody 

c end of cone 

1 index in g ' difectlbn 

j Index in b direction 

o outer boundary 
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t stagnation conditions ■ 

9 >--- f ree-stream conditions . . 

Superscripts: 

n array of known flow properties 
n+1 advanced array of flow properties 

METHOD OF ANALYSIS 


General Equations 

The technique used In the present analysis is a direct adaptation of the 
shock-capturing technique presented In references 4 and. 5. The major details 
of the analysis will be presented here, primarily for reasons of completeness. 
The equations of motion can be written In conservative form,, using vector 
notation and cylindrical coordinates as 

E2 + Fj. + + H = 0 (1) 

where the E, F, G and H vector components are defined as: 
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Equation (1), which represents the continuity and the three momentum equa- 
tions, comprises a complete set when coupled with the energy equation In the 
following form; 

p = p(l - q2) , (2) 

where 

q = /u^ + v^ + w2 (3) 

In applying shock-capturing techniques to the computation of Internal 
flows. It Is most convenient to transform the. physical coordinate system to 
one wherein the distance between the Inner and outer computational boundaries 
Is normalized. This Is done by the following coordinate transformation: 

r - r r - r, . 

Z = z ; 5 = - y ^ ♦ (4) 

o b 
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where, for the case considered here and as shown in figure 1(a), the inner and 
outer boundaries are given respectively by 

Tb = rb<Z) ; to = ro(Z) ; ,(5) 

It should be noted that, for an axially symmetric configuration, the inner 
and outer boundaries are independent of the meridional angle. 

The conservative form of the equations of motion, equation (1), can be 
retained in the transformed coordinate system as: 

: E_ + G_ + +, H = 0 (6) 

^ ^ 

where new variables are defined as follows; ' . 



, A plane of known flow quantities at some constant z” is required to 
initiate the solution. Within this plane, an array of discrete points Is...... 

identified (fig. 1(b)) such that along each meridional ray, there are Nr'; 
points Indexed by 1 < i < N^, and there are meridional rays Indexed by 

1 j N(|). A step size consistent with the stability requirements; of,’ the 
differential equations is first determined for advancing the solution to 
Z^+l. Since only supersonic flows are considered in the present analysis, 
the Z-axis Is the hyperbolic coordinate. For each point in the flow field 
a minimum of four step sizes, corresponding to the four characteristic 
directions, must be found as follows; 



and (9) 

'■ • ■ az = ' ■ - ■ ■ , 

. -■ " 5 ^ • ^ ■ ■■ ; ■ 

The maximum eigenvalues are given by » 
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where , 


and 


V^ '^°(P)]ibcal max; = (l°l 1.^2 D 

‘^q “ ^''^“^•^local max = >»ax (|<TiN|Jp^N|) 


M _ 3g , uv ± a/u^ + - a^ 

1,2 az Ui(u2 - a2) 


N 1-. /uw + a/u2 ,+ w2 - a2\ 

°1,2 - r y (u2 _ a2) J 

H‘4* 1 

Advancing the solution to Z Is controlled by the smallest step size 
determined for the entire array at z”. 


(10) 


( 11 ) 


New conservative variables, Ej 23 4 » are found at step Z^ using 
MacCormick's second-order predictor-corrector technique. Some modification 
to the indexing for these equations is required, depending upon the region of 
the flow being analyzed, as shown in figure 1(a). These predictor difference 
equations are given by , . 


(a) For conical, external, and internal flows, where 1 < i < 


and 1 < j <, N 


- 


<t!’ 


yn +1 


« E 


Ij 


- — fp" 

AC -N 


1 + 1 , j 


- - 


A(|) 




•(b) For conical and external flows, where i =N^ and 1 ^ N^, the 

predicted values are set equal to the free-stream values, which will be 
defined later. 


’ (c)'' For internal flows, where ' 1' = and 1 ^ j ^ N^, 

■ . I ^ . 


where the superbar denotes predicted values at z” . The computations are 
carried out on only one side of the plane of symmetry, 0 < ^ < ti. For the 
predictor step, information at (|i = ir + Aij) is heeded, while (as will be 
seen shortly) information at t = - A(|> is needed for the corrector step. 

This information is obtained by simply mirroring w (l.e., changing the sign 
of w) across the plane of symmetry at the appropriate points and using the 
values of p, p, u, and v from the mirror image point. 


The corrector difference equations are given by 
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(a) For conical and external flows, where 2 < i < Ng and 1 £ j < 


and for internal flows, where 2 i 1 i - 1 and 1 :S j i N^» 

e! 




1 < j < N^. 


7 “ 

Iri-l . AZ / 
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(14) 

conical, external 

, and 

Internal flows 

, where 1=1 and 
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(c) For. Internal flows, where 1 = and 1 < j < N^» 

^ - n+l ^ i f n . M \ 

' * ,4 ■ ' , ' 


A4> 


AC 


(16) 


After both the predictor and the corrector steps, the conseirvatlve 
variables 2 3 4 must be decoded- to. yield: the physical variables. p,p,u,v, 
and w siich tfta£Aew predicted or corrected conservative variables F,G, and H 
can be found. Decoding Is accomplished by the following equations: 


E, 


E. 


. E 2 + ./E2,2 - 4kEj2(l - fc) (1 - v2 - w^) 

= 2Ej(l - k) 


p » — 

u 


1(1 _ u2 - v2 - w2) 


(17) 
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i^ter the advanced array is decoded, the flow along the, centerbody and 
annuliis boundaries will not necessarily satisfy the tangency requirement. As 
described in reference 8, the velocity vector at the wall must be rotated 
through a small angle, Av, and the wall pressure corrected using the . Pfandtl- ’ 
Meyer relationship. The turning angle Av for all j is found as follows: 


(a) For the centerbody. 


Av = sin”^ 


-u. 


.n+1 


tan 6 


n+1 


, /.n+1 


*n+l /. 2 j. 

qj ^tan^ 6^ + Ij 


1/2 


(18) 


where the carets over the symbols denote the decoded values from the corrector 
equations (i.e., those values that do not satisfy tangency). 


(b) For the annulus. 


Av = sin"^ 


/.n+1 »n+l /.n+i 

Ujj tan 6^ -v^^ 






/.n+1 /. 2 X 


n+1 


+ 1 


il/2 


(19) 


The pressure ratio associated with this turning angle is found from 
n+1 


^n+1 

P 


n+1 

= 1 -ypi hv 




Av^ 


( 20 ) 


where : 




Jl 


and 



Note that if Av > 0, an expansion will result, while, if ..Av .< 0, a compres- 
sion will result.’ The new corrected wall pressure, consistent with the tan- 
gency condition, is ’ 





n+1 

2 

/.nfl 

P 


( 21 ) 


and the corresponding density is 



r) 


1 /y j 


This relationship for density assumes that the entropy on the centerbody 
and annulus is everywhere equal to the free-stream values. Once the pressure 
and density are known, the total velocity is found from / 




and the individual velocity components are 


^n+l|-^n+l ■(. g-j 


n+1 .-'ri+l 


(v°^^ - M/n2) 
L 


n+1 >^n+l 
j w 


The quantities L', M, and N are 


= 1 


tan 6 ) 

n2 ) 




« Aii+1 . - , ^n+r 

M ,= u tan 0 +, V 


/tan^ <S + 1 ' . . . j 

equations (24) and (25) can be either the surface :angle of 


The- angle 6 in equations (24) and (25j 
the centerbody or the annulus at 


'• Method of Solution^ , ; . ; , : mv' 

Input to the program consists of the free-stream Mach number, the number 
of points and N,j, (both of which are held constant throughout ? the solu- 

tion), body coordinates (which herein are only a function of ■ Z):, and the 
Ihitlal angle ■ of the conical tip. The remaining free-stream input quantities 
are calculated as follows: 


P. - (v + T 
= (i + T 



(26c) 


(26, d) 
(26e) 
(26f) 

The body coordinates are entered at discrete Z locations. For any 
three points on either the centerbody or annulus, a quadratic of the form 
r = aZ^ + bZ + c was determined. Choice of which coordinates to use for 
defining the coefficients was determined such that for any Z’^'*'^, determined 
from = Z*^ + Az, two of the three points should have Z > . No 

explicit attempt at smoothing the angles at the Input points was made; however, 
the present technique has been found to provide a satisfactory representation 
of 'the actual contours. The slope of any Z location was found by differen- 
tiating the above. 

Solution of the conical flow is obtained by an iterative method. First, 
an outer computational boundary for both the conical and external f lowfleld 
must be found that will allow the shockwave from the cone tip to be captured 
between the centerbody surface and the outer boundary. As shown in figure 1,' 
this boundary is taken to be inclined at an angle, relative to the axis of 
symmetry of the inlet, that is, equal to the angle of attack plus, the planar 
shockwave angle for a flow deflection equal to the cone halfr-angle. Similar 
to the technique of reference 9, the free-stream quantities are assigned 
between the outer boundary and the cone surface at some initial Z start. 

The solution is advanced steps in Z. Since the flow is conical, these 

values can be inserted at Z start and the solution again advanced N^ steps 
in Z, This cycling continues until the pressures on the cone at the end of 
the cycle agree to within an e of 0.0001 of those on the N-1^” cycle, 

the solution was then taken to have converged, and these final flow quantities 
were assigned to the end of the conical flow zone Z^. (The location of Zg . 
is somewhat arbitrary, here taken as Zj, = 1 (see fig. 1(a)), but it should 
always be less than or equal to the portion of the centerbody that can be 
represented as a cone.) Usually, convergence was obtained within eight cycles 
for the cases considered here. 

After the conical flow was defined, the solution was advanced to the 
plane of the annulus lip. Here, the mesh is renoded to provide points 

between the centerbody and the annulus lip. For, the remainder of the solution, 
the outer boundary rQ is the radius of the annulus. ^If the external flow 
outer boundary fell outside the annulus lip, the flow properties at the 
renoded radial points were found using linear interpolation in the external 
flow solution. If the external flow outer boundary fell inside the annulus 
lip, linear interpolation was still used, but free-stream flow quantities 
were assigned between. the external flow outer boundary and the annulus lip. 
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Computation of the internal flow continues until the end of the inputy 
coordinates. No attempt to define the external flow beyond the plane of the 
annulus lip was made in this study. 


RESULTS AND DISCUSSION 


Coordinates of the M = 3.5 axisymmetric supersonic inlet, for which 
shock-capturing solutions Save been obtained herein, are given in table 1, aiid‘ 
shown schematically in figure 2. This inlet, which resulted from the design . 
study of reference 2, has the feature that the minimum area remains in nearly’ 
the same location, relative to the annulus, as the centerbody is translated 
forward for off-design (M^ < 3.5) operation. To obtain maximum performance 
with this inlet, it is necessary to remove the boundary layer in regions of 
strong adverse pressure gradients, such as produced by shockwave boundary- 
layer interactions. Regions of boundary-layer removal remain fixed on the 
annulus but, through a complicated porting system, move rearward on the 
centerbody as it is translated forward. (See ref. 2 for details. ) Since, in 
this analysis, boundary-layer effects and mass exchange through the boundaries 
were not considered, discussion of this boundary-layer removal system is 
important only insofar as indicating regions wherein good agreement between 
theoretical calculations and experimental data should not he expected. ’ : 

Any new computational technique should be' compared with a standard solu- 
tion before being applied to more generalized problems. As shown in refer- 
ence 7, excellent agreement between a discrete-shock, finite-difference tech- 
nique and the method of characteristics was obtained for, planar internal 
flows, with good qualitative agreement also being obtained, with a straight ' 
finite-difference method (shock-capturing technique)^ It remains to be shown 
that similar agreement between a shock-capturing solution and the method of 
characteristics can be obtained for axisymmetric flows, particularly when 
using a shock-capturing, technique developed for three-dimensional flows.. A 
shockwave pattern obtained from the method of characteristics for the inlet at 
= 3.5, a = 0", with the centerbody in the design position (herein defined 
as the position wherein the shockwave from the centerbody tip hits the annulus 
lip at M„ = 3.5) is shown in figure 3, Corresponding static pressure distri- 
butions on the centerbody and annulus, obtained using the method of character- 
istics, are shown in figure 4, along with the pressure distributions obtained-, 
by using the shock-capturing technique. For all shock-capturing solutions 
obtained herein, 20 Intervals in the 5 direction and 10 intervals in the 
direction were used. The shock-capturing solution agrees very well with the 
method of characteristics up to the second shockwave reflection on the center- 
body. This is a fairly^ strong shockwave reflection, with the method of 
characteristics indicating a downstream Mach number less than 1.2. Failure of 
the shock-capturing solution at this point is due to the Mach number becoming 
subsonic in the, Interaction region. Overshoots in static pressure, and hence 
undershoots in Mach number, are typical of the solution of reflecting oblique 
shockwaves using MacCormack's second-order-accurate differencing in a shock- 
capturing technique (ref. 10). Several possibilities exist that can possibly 
resolve this situation: (1) Incorporation of higher-order differencing terms, 

as described in reference 10; (2) incorporation of a damping term near the 
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shockwaves, as was used in reference 11; (3) incorporation of a technique for 
allowing the wall entropy to change across reflection regions; or (4) incor- 
poration of discrete shockwaves into the solution i At present, all four of 
these possibilities are being examined. However, it can.be said at this point 
that until the failure occurs, the solution obtained by the shock-capturing 
technique agrees well with that obtained by the method of characteristics. 


A comparison of the shock-capturing solution discussed above with experi- 
mental dsta is shown in figure 5. Good overall agreement with the experimental 
data. is obtained at a = 0®, except in the regions of boundary-layer removal. 

The experimental data indicate shockwave intersections with the annulus and 
centerbody slightly ahead of the predictions of the shock- capturing technique.; 
This is reasonable, since the presence of boundary-layer effects in the experi- 
mental data would tend to displace the shockwave interactions forward. 

An angle of attack solution (a = 5® for the same centerbody location and 
as the above calculations) is shown in figure 6. Again, the solution 
falls (at station Z/r^ = 4.5) due to the Mach number becoming subsonic in the 
strong Interaction region; nevertheless, some qualitative conclusions concern- 
ing angle-of-attack effects can be drawn. Static pressures along the windward, 

<l> = 0, and leeward;, <|) = tt meridians for both the annulus and centerbody are 
shown. Comparison with the solution shown in figure 5 shows that the magni- , 
tude of the pressure rise across the shockwave reflections is markedly increased 
particularly on the centerbody. Further, evidence that the solution is produc- 
ing a plausible reproduction of the physical flow^^is indicated by the behavior 
of the pressure distribution near the annulus lip on the windward meridian. 

Here, the shockwave from the centerbody tip should pass Inside of the annulus' 
lip and, since the initial angle of the annulus lip is about one degree, a 
local expansion should be generated at the annulus lip. The decrease of pres- 
sure indicated by the solution on the windward meridian near the annulus lip . 
2.86 < Z/rg < 3.1 is in agreement with the postulated physical flow. 

All of the angle-6f-attack experimental data that will be discussed herein 
is for angles of attack of 3®. In the following paragraphs, two different 
Mach numbers, and hence, two sets of different centerbody location data, will 
be discussed, The pattern will be first to show an a = 0® solution for the 
inlet and then an angle of attack solution and corresponding comparison with 
experimental data. Angle-of-attack effects on the internal flow pattern will 
thereby be identified. 

The first ^ch number data that will be discussed is for M^ = 3.3, where 
a forward centerbody translation AZ/r^ = 0.356 was required. Note that even 
a small change in angle of attack requires a change in centerbody position for 
optimum Inlet operation. At this M^ and a = 0®, complete solutions were 
obtained to Z/r^ =5.0, which is beyond the minimum area of the inlet (fig. 7). 
The correspoh4ing solution for a = 3® is shown in figure 8. Although the 
angle-of-attack effects here are not as large as shown in figure 6 , there is an 
indication that the strength of the shock reflections on the leeward side are 
Increased. Here again, a complete solution was not obtained because the Mach 
number became subsonic in the strong reflection region near Z/r^ =4.8 on 
the centerbody. Most of the experimental data are for the leeward side, with 
good qualitative agreement being obtained (except for the viscous effects 
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noted earlier). Some data on the windward side of the annulus were obtained 

and these show the same qualitative behavior' with respect to the leeward data 
as indicated by the shock-capturing solutions. 

The last Mach number data that will. be. discussed is for ^ = • 21 5 i nhere " 
a forward translation AZ/r^ = 0.855 was required. The o' =6“ solution is 
shown in figure 9, where a complete solution beyond the minimum area was 
obtained. At this Mach number, a complete solution was also obtained at angle 
of attack, and again the strength of the shockwave reflections was increased 
(fig. 10). For this case, the centerbody was extended sufficiently far for- 
ward that some expansion' to pressures ' lower than the conical values ‘Vas 
expected, and this is shown by both the' data and the shock-capturing solution. 
Again inyoking the qualifications on the agreement due, to viscous effects, 
good qualitative agreement between the shock-capturing solution and e3q>efimeh^ ' 
tal data was obtained. • ■ 


CONCLUDING RIMARKS 


For all of the, internal flows examined, herein, several observations can ^ 
be made; ^ , .... .... . ... , ....... .. 

,(i) A solution given by the shock-capturing technique agrees well with 
one given by the method of characteristics for a = 0“. ' ' - 

(2) . . becondTorder-accurate differencing can' fail in regions nehr^thd^^^ 
reflection of strong shoclcwaves where the downstream Mach number is hear i.- 
Several possibilities such as higher-order differencing, damping terms hear , 
the shockwaves ,. wall entropy correction, and discrete shockwaves are being 
Investigated to determine if this failure can be overcome. ' ' ’ 


/ ( 3 ) Shbck-^capturlng solutions agree well with experimental data,, except 
in regions of strong viscous effects and boundary-layer removal'. ■. 
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TABLE 1.- Moo = 3.5 INLET CONTOUR DEFINITION 


mm 

r/xa 


r/ra 

Centerbody 

Annulus 

0.0 

0.0 

. ^ 2.86 

1.0 

4.0 

0.70532 

3.1 

1.004188 

4.1 

0.7228 

3.2 

1.0054 • 

4.2 

0.7387 

3.4 

1.0051 

4.3 

0.7512 

3.6 

0.99996 

4.4 • 

0.759 

3.8 

0.9882 

4.5 

0.7625 

4.0 

0.9681 

4.55 

0.763 

4.1 

0.954 

4.6 

0.7625 

4.2 

0.9364 

4.65 

0.7611 

4.25 

0.9261 

4.7 

0.7585 

4.3 ■: 

■ 0.9154 ■ 

4.8 

0.7504 

4.4 

0 . 8949 " 

4.9 

0.7391 

4.5 

0.8768 

5.1 

0.7120 

4.55 

0.8695 

5.3 

0.6829 

4.6 

■ 0.864 

5.5 

0.6525 

4.65 

0.86 

‘ 5.6 ■ 

0.6362 

•• 4.7 • 

0.8572 

5 . 7 '' 

o ; 6 i 8 

4.8 

0.8533 

5.8 

■ 0.5973 

4.9 

0.8511 

5.9 

0.5744 

■ 5.0 

0.8502 

6.0 

, 0.5467 

5.1 

0.85 



. 5 . 6 ’ 

. 0.85 . 



, 5.8 

0.8574 



5.9 

0.8646 



6-0 

0.8735 . 

" . . 

4 * ■ . 
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Figure 4.— Comparison of pressure distributions given by the shock-capturing 
technique and the method of characteristics for = 3.5, a = 0“, and 
design centerbody location, = 0. 




Figure 5.— Comparison of pressure distribution given by the shock-capturing 
technique with experimental ,fpr. vM„ = 3.5, a = 0*, and design centerbody 
location, AZ/rg =0. 
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Figure 10.— Comparison of pressure distributions given by shock-capturing 
technique with experimental data for = 2.5, a = 3®, and 
AZ/r^ = 0.855. 


